Introduction
General remarks
The R pacakge enerscape computes the energy landscape based on models predicting the energy costs of transport. Two models are currently implemented, the ARC model for terrestrial, legged animals from (Pontzer 2016) and the model for cycling expenditures from (Prampero et al. 1979). Adding other models is relatively easy and can be done by writing a custom function in the enerscape_internals.R script. Currently, there is no implementation for user-specified models, but I plan to add this in the future.
An extensive derivation of the ARC model can be found in Pontzer (2016). Here, it is sufficient to say that the ARC model accounts for the energy spent for cross-bridge cycling, generating tension, as well as associated activation-relaxation processes that restore cellular ion gradients - mostly Ca\(^{2+}\) necessary for cross-bridge cycling. Notably, predictions of the ARC model fit very well observed data. In enerscape, the dataset from Pontzer (2016) can be access as the variable pontzer.
library(tidyverse)
library(raster)
library(sf)
library(enerscape)
#' Function to plot the areas
#'
#' @param x is a raster
#' @param poly is a shapefile polygon
plot_area <- function(x, poly, mask = NULL, col = NULL, void = FALSE) {
if (is.null(col)) {
plot(x, axes = !void, box = !void)
} else {
plot(x, col = col, axes = !void, box = !void)
}
if (!is.null(mask)) {
plot(mask, add = T, legend = FALSE,
col = adjustcolor("white", alpha.f = 0.5),
axes = !void,
box = !void)
}
plot(poly$geometry, add = TRUE)
}
#' ARC model predictions
#'
#' @param mass body mass of the animal (kg)
#' @param slope incline of the terrain (degree)
ARC <- function(mass, slope) {
E_ar <- 8 * mass^(-0.34)
E_mec <- 50 * (1 + sin((2 * slope - 74) / 180 * pi)) * mass^(-0.12)
E <- E_ar + E_mec
return(E)
}
pontzer %>%
as_tibble %>%
mutate(ARC = ARC(Mass, Incline)) %>%
ggplot() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
geom_point(aes(ARC, Cost.of.Transport, col = Incline)) +
scale_color_gradient(low = "steelblue", high = "tomato") +
scale_x_log10(n.breaks = 11) +
scale_y_log10(n.breaks = 11) +
xlab("ARC predicitons") +
ylab(expression("E"["COT"]*" (J kg"^"-1"*"m"^"-1"*")")) +
theme_bw()
Computing energy landscapes
To compute energy landscapes using enerscape, we need an elevation rasterLayer, normally called a digital elevation model (DEM). The raster should have equal resolution on both horizontal and vertical axes and a planar coordinate reference systems in meters. The following is a raster that have these characteristics.
## [1] 100 100## CRS arguments:
## +proj=utm +zone=32 +datum=WGS84 +units=m +no_defsWe can then compute the energy landscape of the area for a general animal of 50 kg using the enerscape() function.
## - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)## | Calculating slope## | Calculating work## | Calculating conductance (1 / work)## - Do not use slope with negative values for distance calculationsThe output is an enerscape object, which is a list with elements:
- The neighbors used to compute the transition matrix among cells, either 4 (in chess, rook’s move of 1), 8 (king’s move), or 16 (king’s + knigth’s move).
- The body mass of the animal (kg).
- A rasterStack with DEM, Slope, Work (cost of travel), and Conductance (1 / Work).
- A conductance transition matrix that summarize the conductance of moving between neighboring cells.
The function en_extrapolation() can be used to highlight where ARC predictions extrapolate from the testing data. A rasterLayer masking extrapolations for the incline is returned if available.
## $`Mass extrapolated`
## [1] FALSE
##
## $`Slope extrapolated`
## [1] TRUE
##
## $`Slope extrapolation`
## class : RasterLayer
## dimensions : 367, 509, 186803 (nrow, ncol, ncell)
## resolution : 100, 100 (x, y)
## extent : 852140, 903040, 4659980, 4696680 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
## source : memory
## names : Slope
## values : 1, 1 (min, max)plot(dem)
plot(extr$`Slope extrapolation`,
add = TRUE,
col = adjustcolor("blue", alpha.f = 0.3),
legend = FALSE)
arrows(888523, 4674275, 878791, 4671549)
Least-cost paths and random walks
The algorithm to compute least-cost paths from the R package gdistance is implemented in the en_lcp() function (Etten 2017). en_lcp() takes as input two points for which the least-cost path need to be evaluated.1 Alternatively, random points are generated if simulate_random_points = TRUE:
## 1 of 10 random points ...## 2 of 10 random points ...## 3 of 10 random points ...## 4 of 10 random points ...## 5 of 10 random points ...## 6 of 10 random points ...## 7 of 10 random points ...## 8 of 10 random points ...## 9 of 10 random points ...## 10 of 10 random points ...
The random walk algorithm from gdistance is also implemented as en_passage(). As this may be quite slow, I implemented another way of computing the overall connectivity between two points - or for the whole landscape - using Circuitscape and Omniscape. The initialization files for these to algorithms can be written using circuitscape_skeleton() and omniscape_skeleton(), which create circuitscape.ini and omniscape.ini to be launched from within a Julia terminal. See the example [Landscape connectivity for the Marsican bear (Ursus arctos marsicanus) in the Sirente-Velino Regional Park: Circuitscape and Omniscape implementations.] for details.
Landscape connectivity for the Marsican bear (Ursus arctos marsicanus) in the Sirente-Velino Regional Park.
Here, I assess the overall landscape connectivity for the Mariscan brown bear (Ursus arctos marsicanus) in the Sirente-Velino Regional Park (SVRP) in central Italy parcosirentevelino.it. In particular, I derive the energy landscape for the area and use it as a “resistance matrix” for animal movement, i.e. assuming that the animal would move across low energy landscape trajectories in order to minimize energy costs.
The required input for enerscape using the ARC model are:
- a digital elevation model (DEM) of the area of interest (m)
- body mass of the animal (kg)
I load the DEM and the park polygon data.
## Reading layer `SVRP-polygon' from data source `/home/eb97ziwi/Proj/enerscape-paper/data/orso/SVRP-polygon.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 30 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 857140.5 ymin: 4665004 xmax: 897975.7 ymax: 4691681
## CRS: unknownm <- mask(dem, park, inverse = TRUE)
plot_area(dem, park, m, void = TRUE)
title("Elevation (m)")
scalebar(10000,
type = "bar",
divs = 2,
label = as.character(c(0, 5, 10)),
below = "kilometers")
I calculate the energy landscape for a bear of 140 kg.
## - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)## | Calculating slope## | Calculating work## | Calculating conductance (1 / work)## - Do not use slope with negative values for distance calculations
Or plotting only the cost of travel per cell:
plot_area(en$rasters$Work, park, m, col = topo.colors(100), void = TRUE)
title("Energy cost of travel (kcal)")
I choose two points for which to compute the least-cost path and the connectivity using Circuitscape.
p <- data.frame(x = c(877367, 882653),
y = c(4674192, 4677413))
plot_area(dem, park, m, void = TRUE)
points(p, pch = 20)
And I run Circuitscape in Julia.
julia> using Circuitscape
julia> compute("/home/eb97ziwi/Proj/enerscape-paper/data/circuitscape.ini")## [ Info: 2020-12-19 08:31:24 : Precision used: Double
## [ Info: 2020-12-19 08:31:27 : Reading maps
## [ Info: 2020-12-19 08:31:28 : Resistance/Conductance map has 186803 nodes
## [ Info: 2020-12-19 08:31:36 : Solver used: AMG accelerated by CG
## [ Info: 2020-12-19 08:31:36 : Graph has 186803 nodes, 2 focal points and 1 connected components
## [ Info: 2020-12-19 08:31:36 : Total number of pair solves = 1
## [ Info: 2020-12-19 08:31:37 : Time taken to construct preconditioner = 1.683233535 seconds
## [ Info: 2020-12-19 08:31:37 : Time taken to construct local nodemap = 0.008903604 seconds
## [ Info: 2020-12-19 08:31:37 : Solving pair 1 of 1
## [ Info: 2020-12-19 08:31:38 : Time taken to solve linear system = 0.526888286 seconds
## [ Info: 2020-12-19 08:31:39 : Time taken to write current maps = 0.559248897 seconds
## [ Info: 2020-12-19 08:31:39 : Time taken to complete job = 15.111858446
## 3×3 Array{Float64,2}:
## 0.0 1.0 2.0
## 1.0 0.0 31.8625
## 2.0 31.8625 0.0Sometimes, few cells with extremely high current/connectivity can make the output map hard to read, as most of the connectivity values are squeezed into a narrow range. To avoid this, the value of extremely high outliers can be forced to an upper limit. For example, I force outliers to the 0.999 quantile value of the connectivity and standardize the Circuitscape output between 0 and 1.
cs <- raster("../data/orso/circuitscape_cum_curmap.asc")
cs <- log10(cs)
cs <- cs - min(values(cs), na.rm = TRUE)
cs[cs >= quantile(cs, 0.999)] <- quantile(cs, 0.999) #remove far outliers
cs <- cs / max(values(cs), na.rm = TRUE)
plot_area(cs,
park,
m,
col = c("black", viridis::viridis(99)),
void = TRUE)
points(p, pch = 20, col = "black")
contour(mask(dem, m, inverse = TRUE),
add = TRUE,
nlevels = 5,
lt = 2,
lw = 1.5,
col = "grey80")
lines(lcp$Path)
title("Connectivity between locations")
To compute the overall landscape connectivity, I use the Omniscape algorithm, which runs Circuitscape iteratively on a moving window. I initialize the omniscape.ini file to be passed to Julia.
omniscape_skeleton(en,
file.path("/home",
"eb97ziwi",
"Proj",
"enerscape-paper",
"data"),
radius = 10)And I run the Omniscape algorithm in a Julia terminal. Don’t forget to specify the number of parallel threads: https://docs.circuitscape.org/Omniscape.jl/stable/usage/#Parallel-Processing.
julia> using Omniscape
julia> run_omniscape("/home/eb97ziwi/Proj/enerscape-paper/data/omniscape.ini")## Starting up Omniscape. Using 7 workers in parallel. Using double precision...
## Solving moving window targets...
## Progress: 100% | Time: 0:03:14
## Done!
## Time taken to complete job: 196.3408 seconds
## Outputs written to
## /home/eb97ziwi//home/eb97ziwi/Proj/enerscape-paper/data/omniscape_2
## (Union{Missing, Float64}[0.5152914765395912 0.5884433266779532 … 0.719651357069539 0.8600629588479131;
## 0.4834003129933446 0.6594116484239914 … 0.635221929581483 0.7110678186022065;
## … ;
## 0.2127925084620419 0.2067469779201847 … 0.39577453680538544 0.46681040575465904;
## 0.17382917047741175 0.2093726458198999 … 0.31763823508216504 0.3738227777731409]
## Union{Missing, Float64}[1.0885252650071 1.136985249629082 … 1.169180462521136 1.1104937306261533;
## 1.053217843455199 1.0090282247541849 … 0.9301240063409788 1.1859951520336032;
## … ;
## 1.1798282079383584 0.9871029439774996 … 0.9873244151514602 1.2432973368603717;
## 1.1120398150491893 1.168088759875787 … 1.0583987239714665 1.0877633126993709])I remove outliers and standardize the Omniscape output between 0 and 1.
os <- raster("../data/orso/omniscape/cum_currmap.tif")
os <- os - min(values(os), na.rm = TRUE)
os[os >= quantile(os, 0.999)] <- quantile(os, 0.999) #remove outliers
os <- os / max(values(os), na.rm = TRUE)
plot_area(os,
park,
m,
col = c("black", viridis::cividis(99)),
void = TRUE)
contour(mask(dem, m, inverse = TRUE),
add = TRUE,
nlevels = 5,
lt = 2,
lw = 1.5,
col = "grey80")
title("Landscape connectivity")
Habitat use of feral animals in the “Rewilding Mols” experiment in the Mols Bjerge National Park.
Here, I assess how the energy costs of locomotion affect differently Galloway cattle and Exmoor ponies in a rewilding experiment in Denmark. Energy landscapes are used to compute the landscape connectivity of the rewilding site, which is then used to model species habitat usage. First, I load the DEM and associated layers.
## Reading layer `Mols' from data source `/home/eb97ziwi/Proj/enerscape-paper/data/mols/Mols.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 596616 ymin: 6231887 xmax: 598918.6 ymax: 6233212
## CRS: unknown## Reading layer `Mols' from data source `/home/eb97ziwi/Proj/enerscape-paper/data/mols/Mols.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 596616 ymin: 6231887 xmax: 598918.6 ymax: 6233212
## CRS: unknowndem <- raster("../data/mols/DK-DEM.asc") %>%
crop(extent(st_buffer(park, 100)))
crs(dem) <- crs(park)
m <- mask(dem, park, inverse = TRUE)
plot_area(dem, park, m, void = TRUE)
I compute the energy landscapes for both species and calculate the landscape connectivity using the omniscape algorithm.
en_horse <- enerscape(dem, 320, "kcal")
omniscape_skeleton(en_horse,
"/home/eb97ziwi/Proj/enerscape-paper/data/mols/omniscape_horse/",
radius = 10)
en_cow <- enerscape(dem, 550, "kcal")
omniscape_skeleton(en_cow,
"/home/eb97ziwi/Proj/enerscape-paper/data/mols/omniscape_cattle/",
radius = 10)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition
I use landscape connectivity as an explanatory variable in a RandomForest (RF) model to predict the density of GPS records for both species, including as other explanatory variables NDVI, tree cover density, and the other species GPS record density. I skip here RF models and load directly the results of the analysis.
horse <- read_rds("../data/horse_rf.rds")
cattle <- read_rds("../data/mols/cattle_rf.rds")
# inspect table of RF model
head(horse$data)## # A tibble: 6 x 8
## n ndvi tree omni cattle month lon lat
## <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 -0.0849 0.139 81.7 0.139 0.521 Jan 597395 6233175
## 2 0.0251 0.117 74.5 0.162 0.612 Jan 597425 6233175
## 3 0.171 0.143 75.1 0.272 0.954 Jan 597395 6233145
## 4 0.356 0.148 84.1 0.338 1.09 Jan 597425 6233145
## 5 0.528 0.115 78.9 0.164 1.08 Jan 597455 6233145
## 6 0.429 0.152 64.6 0.399 1.11 Jan 597395 6233115I plot the importance of predictor variables for both species.
importance <- tibble(
var = c(rownames(horse$importance), rownames(cattle$importance)),
IncNodePurity = c(horse$importance, cattle$importance),
species = c(rep("Exmoor pony", 5), rep("Galloway cattle", 5))
) %>%
mutate(var = modify(var, function(x) {
if (grepl("GPS", x)) {
"Density of the other species"
} else {
x
}
}))
importance %>%
ggplot() +
geom_point(aes(reorder(var, IncNodePurity), IncNodePurity, col = species),
position = position_dodge(width = 0.5)) +
coord_flip() +
scale_color_manual(name = "", values = RColorBrewer::brewer.pal(3, "Set2")) +
xlab("") +
theme_classic() +
theme(panel.grid.major.y = element_line(color = adjustcolor("grey20", alpha.f = 0.2),
linetype = "dashed"))
The variable importance plot shows that landscape connectivity (omni) strongly influence cattle habitat usage, but it has only limited effects on the area used by horses. Partial plots show that cattle tend to use habitats that are well-connected, i.e. that can be reached with low energy expenditures.
horse <- read_rds("../data/mols/horse_partialplots.rds")
cattle <- read_rds("../data/mols/cattle_partialplots.rds")
res <- tibble()
for (vars in names(horse$plot)) {
res <- res %>%
bind_rows(mutate(horse$plot[[vars]],
var = vars,
species = "Exmoor pony"))
}
for (vars in names(cattle$plot)) {
res <- res %>%
bind_rows(mutate(cattle$plot[[vars]],
var = vars,
species = "Galloway cattle"))
}
res %>%
mutate(var = modify(var, function(x) {
if (grepl("GPS", x)) {
"Density of the other species"
} else {
x
}
})) %>%
ggplot() +
aes(x, y, col = species) +
geom_line(size = 0.5) +
geom_point(shape = 21, size = 0.75, fill = "white") +
scale_color_manual(name = "", values = RColorBrewer::brewer.pal(3, "Set2")) +
facet_wrap(~var, scales = "free") +
theme_classic() +
xlab("") +
ylab("Marginal effect on GPS KDE") +
theme(strip.background = element_blank(),
strip.text = element_text(hjust = 0.5))
Overall, these findings highlight that energy costs of locomotion influence the habitat usage for cattle, but not for horses, in the rewilding experiment, suggesting the expected ecological effects of species to act differently accoring to topographical features of the site.
Energy costs of cycling L’Eroica
Here, I study how body weight (kg) and cycling speed (km/h) influence the energy required by participants to cycle L’Eroica event. In particular, I derive the energy landscape for the route and sum the energy costs of all cells.
I load the DEM and the route polygon data.
dem <- raster("../data/eroica/eroica.tif")
route <- st_read("../data/eroica/route.shp") %>%
st_transform(crs(dem))## Reading layer `route' from data source `/home/eb97ziwi/Proj/enerscape-paper/data/eroica/route.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 12 fields
## geometry type: LINESTRING
## dimension: XYZ
## bbox: xmin: 11.3464 ymin: 43.03679 xmax: 11.61873 ymax: 43.49132
## z_range: zmin: 129.4 zmax: 620
## geographic CRS: WGS 84r <- st_coordinates(route)[, 1:2] %>%
Line() %>%
Lines("Eroica")
r <- SpatialLines(list(r))
crs(r) <- crs(route)
plot_area(dem, route, void = TRUE)
scalebar(d = 10000,
type = "bar",
below = "kilometers",
label = as.character(c(0, 5000, 10000) / 1000),
div = 4)
I force the DEM on the route only.
The energy landscape for a cyclist of 60 kg cycling at an overall speed of 30 km/h can then be computed using to use the model from (Prampero et al. 1979) by specifying method = "cycling".
## - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)## | Calculating slope## | Calculating work## | Calculating conductance (1 / work)## - Do not use slope with negative values for distance calculationsThe total energy costs of cycling L’Eroica are then derived by summing the energy expenditures for the whole route.
## [1] 5782.049To get the energy costs for cyclists of different weights and for other cycling speed, I make a for loop.
library(foreach)
res <- foreach (weight = seq(50, 100, by = 10), .combine = "rbind") %do% {
res <- foreach (speed = seq(10, 50, by = 5),
.combine = "rbind") %dopar% {
en <- enerscape(el,
weight,
"kcal",
4,
method = "cycling",
v = speed)
cbind(speed,
sum(values(mask(en$rasters$Work, el)),
na.rm = TRUE))
}
cbind(weight, res)
} %>%
as_tibble()I can then plot the energy costs for different body weight and cycling speed.
res %>%
mutate(weight = as.factor(weight)) %>%
ggplot() +
geom_line(aes(speed, V3, col = weight)) +
geom_violin(aes(speed, V3, group = as.factor(speed)), alpha = 0, draw_quantiles = 0.5) +
geom_point(aes(speed, V3, col = weight), size = 1) +
xlab("Speed (km / h)") +
ylab("Energy consumed (kcal)") +
scale_color_brewer("Weight (kg)", palette = "RdGy") +
theme_classic() +
theme(panel.grid.major.y = element_line(color = "grey80", size = 0.25))
and the difference in minimum and maximum costs at same speed, but different weight.
res %>%
group_by(speed) %>%
summarize(Min = min(V3),
Max = max(V3)) %>%
mutate(Diff = Max - Min) %>%
ggplot() +
aes(speed, Diff) +
geom_smooth(alpha = 0.2, col = "tomato") +
xlab("Speed (km / h)") +
ylab("Maximum energy costs\ndifference (kcal)") +
geom_point() +
theme_classic()## `summarise()` ungrouping output (override with `.groups` argument)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
References
Etten, Jacob van. 2017. “R Package Gdistance: Distances and Routes on Geographical Grids.” Journal of Statistical Software 76 (1): 1–21. https://doi.org/10.18637/jss.v076.i13.
Pontzer, Herman. 2016. “A Unified Theory for the Energy Cost of Legged Locomotion.” Biology Letters 12 (2): 20150935. https://doi.org/10.1098/rsbl.2015.0935.
Prampero, P. E. di, G. Cortili, P. Mognoni, and F. Saibene. 1979. “Equation of Motion of a Cyclist.” Journal of Applied Physiology 47 (1): 201–6. https://doi.org/10.1152/jappl.1979.47.1.201.
A path can also be drawn on map using
en_path()specifying draw = TRUE and the number of points to draw the path, e.g.(en, draw = TRUE, n = 10).↩︎